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Abstract. The Escalator Boxcar Train method (EBT) is a numerical method 
for structured population models of McKendrick - von Foerster type. Those 
models consist of a certain class of hyperbolic partial differential equations and 
describe time evolution of the distribution density of the structure variable de¬ 
scribing a feature of individuals in the population. The method was introduced 
in late eighties and widely used in theoretical biology, but its convergence was 
proven only in recent years using the framework of measure-valued solutions. 

Till now the EBT method was developed only for scalar equation models. In 
this paper we derive a full numerical EBT scheme for age-structured, two-sex 
population model (Fredrickson-Hoppensteadt model), which consists of three 
coupled hyperbolic partial differential equations with nonlocal boundary con¬ 
ditions. It is the first step towards extending the EBT method to systems of 
structured population equations. 

1. Introduction. 

1.1. Problem formulation. The Escalator Boxcar Train (EBT) algorithm is a 
numerical method for solving structured population models. It is based on repre¬ 
senting the solution as a sum of masses localised in discrete points and tracing their 
evolution due to transport and growth. The algorithm has been used in applications 
for a long time [4], however convergence of the scheme was shown only recently in 
ref. [1] (without the rate of the convergence), and later in ref. [6] (the rate of the 
convergence) using the metric space approach proposed in ref. [7, 8]. So far the 
method has been established only for single equation models. The aim of this paper 
is to derive the EBT algorithm for a system of structured population equations. 
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We focus on the Fredrickson-Hoppensteadt model, which is a two-sex population 
model with age as a structure variable. The model was originally formulated in 
ref. [5] and later developed in ref. [11]. The model consists of three population 
equations with structure coupled through nonlocal boundary terms and a nonlocal 
and nonlinear source term. Functions x) and {t, x) describe the distribution 
of males and females, while x, y) is the number of couples at time t. Structural 
variables x and y denote the age of males and females, respectively. The following 
system of equations describes the dynamics of the population of males, females and 
couples: 


dtu™'{t, x) + dxu"^{t, x) -\- c™(f, x)u™‘{t, x) 

u’"(0,a:) 



dtu-f{t,y) +dyU-f(t,y) + c-^ (t,y)u^ {t,y) 

0 ) 

{0, x) 


0 , 

/ P^{t,x,y)u^{t,x,y)dxdy, 

ul{x), 

( 1 ) 


X, y) -k X, y) + dyu‘^{t, x, y) -k c^(t, x, y)u°{t, x, y) 

0) = u°(f,0, y) 

u'=(0,a:,y) 


T{t,x,y), 

0 , 

u^Q{x,y). 


Functions c™, and c‘^ describe the rate of disappearance of individuals. In the 
case of males and females disappearance equals death, while in the case of couples 
it equals divorce or death of one of spouses. Functions /3™ and are birth rates 
of males and females. In the general case the functions may depend nonlinearly on 
ecological pressure, which is usually modelled by a nonlocal operator depending on 
the distribution of males, females and couples. However, in this paper we restrict 
our considerations to the case where disappearance and birth rates depend only on 
time and on the structure variable. 

The marriage function 'T models the number of new marriages of males and 
females of age x and y, respectively, at time t. It depends nonlinearly on the distri¬ 
bution of individuals. The choice of this function is a subject of ongoing discussions, 
see [9], [10], [16], [13], due to the properties like heterosexuality, homogeneity, consis¬ 
tency or competition. In this paper we follow the formulation proposed in ref. [12], 
namely: 


T{t,x,y) = F(t,u'^{t,x),uf(t,x),u^{t,x,y)) = (2) 

Q{x, y)h{x)g{y) [u^{t, x) - u‘^{t, x, y)dy] y) - /” x, y)dx\ 

O' + ir x) - X, y)dy] dx + /“ g{y) [uf{t, y) - /“ x, y)dx] dy 

Function Q{x,y) describes the marriage rate of males of age x and females of 
age y. Notice that [u™(f,x) — u'^(t,x,y)dy1 is the amount of unmarried males 
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and [u^ {t,y) — u‘^{t,x,y)dx\ is the number of unmarried females. Function 

h € L^(R+) n L°°(R+) describes the distribution of eligible males on the marriage 
market. Function g is of the same regularity and describes the distribution of eligible 
females on the marriage market. 

1.2. EBT method for a scalar equation. Before moving to the derivation of 
the numerical scheme for system (1), we briefly discuss the analytical and numerical 
results for a simpler model of the McKendrick - von Foerster type [14] : 

dtu{t,x) + dx(b{t,x)u(t,x)) + c{t,x)u{t,x) = 0, (3) 

pOO 

u{t,Xb) = / I3{t,x)u(t,x)dx, 

Jxf, 

u{0,x) = uo{x), 


which is a one-sex population model with structure variable x. The function b{t, x) 
is the rate at which the individuals change their state. In particular case, when 
X represents age, function 6(t, x) equals one. One of the many numerical methods 
applied to this problem is the Escalator Boxcar Train introduced in ref. [4]. 

EBT method is based on representing the solution as a sum of masses (cohorts) 
localised in discrete points and tracing their spatio-temporal evolution using the 
following algorithm: 


d . . 

= 

b{t,x,{t)), 

for i = B + 1,.. 

,J, 

= 

-c{t,Xi{t))mi{t), 

for i = B + 1,.. 


/ 

Xsit) 

( ^B(t) , 

if mB{t) > 0, 


\ Xb, 

otherwise. 


d 

= b(t,Xb)mB{t) + 

dxb{t,Xb)TrB (t) 


-c{t,Xb)TTBit), 



d , , 

= -c{t,Xb)mB{t) 

- dxC(t,Xb)TTB{t) 






= 0, 



71^(4) 

= 0. 




The procedure consists in solving the system of ordinary differential equations 
(ODEs) (4)-(5) on a sufficiently short time interval [tk, tfe-i-i], iteratively with respect 
to k = 0,1,.... The details of the EBT algorithm are discussed in Section 2 for 
system (1) and therefore we omit them at this stage. 

As mentioned before, the rate of convergence of the method for a nonlinear 
version of model (3) has been recently shown in ref. [6]. The results are based on 
Lipschitz semiflows in metric spaces developed in ref. [3] and applied to structured 
population models in ref. [8]. 
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Definition 1 . 1 . Let {E,p) be a metric space. A Lipschitz semiflow is a semigroup 
S : [0, (5] X [0, T] X A —>■ A satisfying 

p{S{t; t)p,S{s-, r)i/) < L{p{p, ly) + \t - s|), 

where s,t G [0, 5], r, s + r, t + t € [0, T] and G E. 

The idea of the proof of convergence presented in ref. [6] consists in considering 
the system (3) (under the assumption that b,c,P G C“([0,T]; W^’°°(R+)) as a 
semiflow in a space of nonnegative Radon measures (Af"^(M+)) equipped with the 
Lipschitz bounded distance, defined below. 

Definition 1 . 2 . Let G A1'^(R.+ ). The distance function pp : A1'''(M+) x 
A1 ■*■(]&+) ^ [0,oo] is defined by 

Pf{p,v) :=swp[ j ipd{p-v)\'ip GC\R+),\\il}\\n'i,,^ KlV ( 6 ) 

Jr+ '' 

where ||i/)||v^i.=o = max{||i/;||oo, ||9a;'0||oo}- 

The metric pp is a distance derived from the dual norm of with the space 

pj/i,oo equipped with its usual norm, i.e. 

||7||vvi,=o = max{||7||L=o, \\d^j\\L--}, 

see ref. [15] and [18]. 

The solution to the system (4)-(5) is interpreted as an element of the space of 
nonnegative Radon measures given by the formula 

.7 

B(t) 

which means that /j," (f) is a linear combination of Dirac deltas with masses localised 
in points Xi{t). 

The proof of convergence of the ETB scheme is based on the following proposi¬ 
tion, the proof of which is a modihcation of the proof presented in ref. [2] . 

Proposition 1 . Let S : E x [0,(5] x [0,T] -G E be a Lipschitz semiflow. For every 
Lipschitz continuous map [0, T\ ^ t ^ vt G {E, p) the following estimate holds: 

p{iyt,S{t;0)po) < L f liminf dr, (8) 

J[o,T] h 

where p is a corresponding metric. 

Application of Proposition 1 reduces the proof of the convergence of the measure 
given by formula (7) to the solution of the model (3) to the estimate of the right 
hand side of inequality (8). In ref. [6] it was proven analytically and confirmed 
numerically that the method is of the first order. 

1.3. Organisation of the paper. This paper is devoted to derivation of a corre¬ 
sponding EBT method for the two-sex age-structured population model (1). The 
convergence of the method follows from stability results obtained in ref. [17] and its 
proof is deferred to a forthcoming paper. Section 2 is devoted to derivation of the 
EBT algorithm for the two-sex population model. In Section 3, we provide well- 
posedness results for the structured population model and for the EBT scheme, and 
discuss a sketch of the convergence proof. 
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2. Derivation of the EBT algorithm for the two-sex population model. 

This section is devoted to the formal derivation of the EBT scheme for system (1). 
We assume that the solution u-^, i u'^ to system (1) are continuously differentiable 
with respect to time and twice time continuously differentiable with respect to 
the structure variables. Additionally, the solution itself and its first derivative 
with respect to time and first and second derivatives with respect to the structure 
variables are assumed to be bounded. The same regularity is assumed on all model 
coefficients appearing in the underlying system (1). 

The first step of the EBT method is grouping the initial distribution of males, 
females and couples into the so-called cohorts, with the grouping performed with 
respect to the age of the individuals. While in the case of males and females the 
groups form one-dimensional intervals including individuals of certain age, in the 
case of couples we deal with two-dimensional intervals to take into account the age 
of both partners. Each cohort is characterised by two quantities: the mass and its 
location. The mass (m'"(t), is the number of individuals or cou¬ 

ples within the cohort at time t, and its location (a;"*(t), y^{t) and (a:°(t), 
respectively) is the average value of the structure variable over the underlying co¬ 
hort. Cohorts are evolving in time along the characteristic lines of the underlying 
problem. We wish to track how masses and locations of cohorts change in time. 
Once individuals or couples are assigned to a certain cohort they remain there till 
disappearance (death or split-up in case of couples). The influx of new couples is de¬ 
scribed with marriage function, while the influx of new males and females is given 
in the boundary conditions. Each time step can be treated as an internalisation 
moment, where a new boundary cohort appears. In case of couples the boundary 
cohorts are empty and in case of individuals the masses in boundary cohorts result 
from the boundary conditions. 

Following this argument, we start with introducing the initial cohorts. At time 
t = 0 we impose space mesh points ^("(0) and (0), where i, j = Bq, ..., J in such 

a way that = IbqW ~ define J — Bq initial intervals for males and 

females 


f^m(0) = r(0),/.^i(0)). f^^+i(0) = 






hj = Bo,..., J - 1, 


and [J — Bq)^ two-dimensional intervals for couples 

fl(,+i)(,+i)(0) = 0-1 (0) X Oj+i(0), *, J = Bo,..., J - 1. 

The above partitioning requires the condition that the supports of initial functions 

f 

M™ , Uq and which are the initial distributions of individuals and couples, are 
contained in the sums of all corresponding cohorts that is 

J J 

supp«) C O-(0) = U O-(0), supp(u^) C 0^(0) = y O^^(O), 

i—S q + I j=BQ-\-l 


supp«) C 0^(0) = y O“^(0). 

*.i=Bo-|-l 

Bearing in mind that the cohorts evolve in time, we write {Oi(t)}/,^^^_|_i, 

and understand that functions and l{ (t) are 
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characteristics defined by the transport operators, thus they satisfy differential equa¬ 
tions = 1 and = 1 and are straight lines: 


lT{t) = t + lT{0), i = 
ijit) = t + Z/(0), 


We impose a mesh on the time variable t € [0,T) in such a way that to = 0 and 
^"+i) = 1^0'^ 0 < t < fi, we define boundary cohorts = 

[OJBoit)), = [Oji^it)) and X = ^Ti^) X 

X where i,j = Bq + 1, J. Mesh points where 

tn = 1,..., iV, are called internalisation moments. 

Altogether, we have J — Bq + 1 initial cohorts for males and females and (J — 
Bq + I)'^ initial cohorts of couples. 

At each point tn, we define Bn = Bq — n and (t) = t — tn for tn < t < t„+i. 
Now, for t S [tn,tn+i), we define internal cohorts: 


o: 

n 


i+l(^) 
/ 


and boundary cohorts: 


forz = S„_i...,J-l 
for 

^T+i{t) X nf+,{t), for 1 ,J = Bn-i ..., J - 1 


^bA*) = AIbA*)) 

= n^At)x^iAi) 

(t) = Or(t) X (t), for i = Bn-i ..., J -1 

it) = ^B^ it) X nt {t),fOTJ=Bn-l...,J-l. 


This means that at each internalisation moment we are adding new boundary 
cohorts to account for the influx of new individuals. At the same time boundary 
cohorts from the previous time interval [tn-i,tn) became internal ones on 

It is essential to observe that at time step we deal with J — Bq + n + 1 male 
and female cohorts and {J — Bq + n + 1)'^ couple cohorts. 

For tn < t < tn+i masses and their locations are defined as follows: 


r Birw = 
^sAt) = 
[ ^Tit) = 
toJ (t) = 

' viAt) = 

Vj it) = 

V 


x)dx, i = Bn, ..., J 
f 0 if it) = 0, 

1 J otherwise 

mfTty ^U^it^ x)dx, i = Bn + I, ■ ■ ■ , J 

'^^it’ y)^y^ j = Bn,..., J 
r 0 if mB^{t) = 0, 

(t) " , 

—, otherwise 

I ™k(‘) 

/n-(t) y^^it’ y)dy, J = Bn + l,...,J 


(9) 


(10) 
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liif. it) 2 ;, y)dxdy, i,j=Bn,...,J, 


(11) 


where n^^(t)= f^^^xu^{t,x)dx and {t,y)dy. 

Let us notice that we deal with 3 unknowns, thus 3 equations, for each male and 
female cohort and with 2 unknowns for each couple cohort. Knowing the amount 
of cohorts explicitly we deduce that at each time step we have 3{J — Bq + u+I) 
equations for males and females and 2(J — Bq + n + l)^ equations for couples. 
Consequently, at we have 8 + 2(J — i?o + n) equations more then we 

had at t„_i <t < in¬ 
ks our aim is to track the dynamics of the above functions describing masses and 
their locations, we wish to derive the appropriate system of ordinary differential 
equations. To shorten the notation, from now on, we assume that t G [tmtn-\-i) 
and we write B instead of Bn- Let us start with formulating ODEs for masses and 
locations of male individuals and symmetric results for female individuals and then 
move to the more demanding case of the couples. We start with differentiating 
xY^{t), i = B + I, - - -, J, and boundary functions mg{t) and n^(t): 


dt 


\t) = 


[ dtu^{t,x)dx + ^li{t)u"^{t,li{t)) - 

-iQ-f-it) ® 

f x)dx + hit)) — li-i{t)) 

-JQY'it) 

( dtu™'{t,x)dx + f dxU^{t,x)dx 

— ( c^{t,x)u^{t,x)dx, 

jQ.'^it) 


dt 


m = 


In'S it) 

[ dtu"^(t, x)dx + ^lB{t)u"^{t,lB{,t)) - U^{t, 0) + u^{t, 0) 

Jnsit) dt 

[ dtu"^{t,x)dx + f dxu'^{t,x)dx + f I3"^{t,x,y)u^{t,x,y)dxdy, 
Jnvi(t) Jn’s(t) 


dtu'^{t,x)dx + —lB{t)u'^ {t,lB{t)) 


/ c^{t,x)u"^(t^x)dx + / l3"^{t,x,y)u'^(t^x,y)dxdy- 

/n™(t) Jr^ 


To differentiate the location function a;™(t), for z = i? + 1,..., J, we check how 
the first moment xY^{t) = xu"^(t,x)dx evolves. 
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dt 


xTit) = 


[ xdtu""{t, x)dx + li-i{t)) 

JqY'W 

f xdtu'^{t,x)dx + f dx{xu'^{t,x))dx 

f xdtu'^{t,x)dx + ( xdxU™‘{t,x)dx + f u'^{t,x)dx 

{t,x)u'^{t,x)dx + / u^{t,x)dx. 


a™(t) 


Having calculated, it is easy to differentiate the location functions for 

the internal cohorts: 




cTit) 


j-^xT{t)mTit) xrit)irnr{t) 


dt 

—^ / aic™ (t, x)u"^ (t, x)dx + I u^(t, x)dx 

Hi dn™(t) Hi JQY'W 


dt 

1 


1 1 


m™(t) m™(t) 


{t, x)dx / c™ (t, x)vr^ (t, x)dx 
JQY'it) 


=x^{t) 


i.,', / / u'^(t,x)dx 


We also need to differentiate the first moment, n^(t), in the boundary cohort: 


dt 


U^{t) = 


dt 


Q,b {t) 


xu^(t, x)dx 


xdtu^{t, x)dx + lB{t)u'^{t,lB(t)) — 0 • u^{t, 0) + 0 ■ u"‘(t, 0) 

xdtu’^(t,x)dx + f dx{xu'^{t,x))dx, 

Jng(t) 

/ xdtu'^{t,x)dx + / xdxU^{t,x)dx + / u"^{t,x)dx, 

[ xc^{t,x)u^{t,x)dx + [ u^{t,x)dx. 

Jnvi(t) JQ’S(t) 


'n’git) 

[ 

In'Sit) 

f 

In'Sit) 

L^it) 


Because we need to obtain the closed form of the schemes, we use the following 
facts that hold for i = H + 1,..., J: 
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— x)u"^{t, x)dx 
f{x)u'^{t, x)dx 


0 , ( 12 ) 

+ O f max \x - a;r(i)p)l,3) 


for / G C2(R+). 

While the first one can be easily justified based on the definition of x'^{t) 

- x]u^{t,x)dx = a;rW/n.(t) u"^{t,x)dx - xu^{t,x)dx 
= xY^{t)m^(t) - m™(t)a;™(t) = 0, 

the proof of the second one requires equation (12) and the first order Taylor ap¬ 
proximation 



/(x)u™(t, x)dx 


a™(t) 


f{x"^{t))u"^{t, x)dx 


+ [ [x-xY^it)]u"^(t,x)dx 

u-x JQY'it) 

+ [ O {\x — x'^{t)^) u'^{t,x)dx 
JnT(t) 


f{xT(t))mT{t) + ||u'"(t, OliLqn-p)) ■ O I max |a; - xY^{t)\- 

\aiGr2^(i} 


Furthermore, we observe that 



X-X'r^{t) 

y-vtjit) 



h{x,y) 

f2{x,y) 


u^{t, X, y)dxdy 


x, y)dxdy 


\ h{x%{t).yt,{t)) j 

->rO \ max I {x1 (t) - X, y^, (t) - y) M 


(14) 

(15) 


for/i,/2GC2(R2). 


Equality (15) holds because 
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X, y)dxdy 


% 

dx 

a/2 

dx 


dy 

0/2 

dy 


{xUt),yUt)) 


y-yf-iit) 


X, y)dxdy 


+ I 0{ miix ^{t)-x,y^ {t)-y)\'^ \u‘^it,x,y)dxdy 

\[x,y)&Q.Ut) ■> ] 


/ fiixijit)^ytj(t)) 

V f2{x%{t),ytj{t)) 


mlAt) 


+ Ol max \{Xij(t) - x,y^j{t) - y)\^ ■ •, •)||Li(nf.(t)) 


Now we are ready to derive the system of equations for the male population: 


^rn^{t) = — f c^{t,x)u'^{t,x)dx 

dt Jn^it) 

= -c^{t,xT{t))mT{t)+o( max \x-xT{t)A, ( 16 ) 

/ [x'^{t)-x\c'^{t,x)u^{t,x)dx 

dt TOj (t) JciY^it) 

+ 3 ^ 77 ^ / u"^{t,x)dx ( 17 ) 

= - xT{t)]c"^{t,xT{t))mT{t) 

+° CiSf - i ,, I * - ‘'^<‘>'0 • 
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f c^{t,x)u"^{t,x)dx + [ j3‘^{t,x,y)u^{t,x,y)dxdy 
ing(t) 

f c^{t,0)u'^{t,x)dx — ( dxC^{t,0){x — 0)u'^{t,x)dz 

O [\x\'^) u"^{t,x)dx + ( (3"^{t,x,y)u‘^{t,x,y)dxdy 
jRi 


d) 

= -c™(t,0)m5(i)-a,c™(t,0)ng(<) + C> f max \xf 

Vxeng(t) 


+ ^ r{t,x<r^{t),yt,{t))m^,{t) 


(18) 


iJ—B 

J 


+ XI , max |(a;)l(t)-a;,y))(i)-y)|^ 

i,j=B 


= -C™(t,0)m5(i)-a^c™(t,0)ng(t) + C) max |a;| 

\xGQ^{t) 

i,j=B 

+ max 0\ max {{x^Jt) - x,yfJt) - y)\'^ ] 

r,3=B,...,J I (x,y)en9.(i)" 


and finally 


_d 

dt 


= ~ [ xc'^{t,x)u'^{t,x)dx + f u^{t,x)d. 


Q’Sit) 


= (i) — c'"(f, O)n^(f) + O ( max \x\^ 

\xeQ^{t) 


(19) 


As we mentioned before, one can easily conclude that the system of ODEs for the 
female population is of similar form, namely 


s"'*') = - 


m(t) 


c^{t,y)u^t,y)dy 


= -cHt, Vj {t) + o[ max \y - yj(t)p 

Xyeiiiit) 


(20) 
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mUt) JnUt) 


+ 


mUt) JaUt) 


it, y)dy 


rr^Wj it) - yj it)Wit, yf it) 


mt it) 


+ (^) + ^ ( \y - Vj (^) 1 ^ ) 

mt it) \y^^i W / 


= 1 + 0 \ max \y - yf (i)p ) , 


\yent (t) 


(21) 


d 

dt 


T^sit) = - ct{t,y)ut^it,y)dy+ / pt' it,x,y)u°it,x,y)dxdy 

= - f c-t{t, 0 )u-tit,y)dy - [ it,Q)iy - D)ut it,y)dy 

+ [ O i\y\'^) ut"it,y)dy + [ it,x,y)u‘^it,x,y)dxdy 

Jni(t) JR^. 


= -ctit,0)mtgit) - d^ctit,0)lltgit) + O \ max |j/|^ 


yen(j(t) 


+ l^^it,^ijit),yijit))'mtjit) 


( 22 ) 


i,j=B 

J 


+ ll«°(^^^')llLi(n=^(t))0 max |(a:=,(t)-x,y'j(t)-y)|2 

i,j=B 

= -ctit,0)m-tgit) - da:ct{t,0)ll-t^it)+0 i max |yM 

J 

j 

+ l^^it,^+it),y'rjit))mdljit) 

i,3=B 

+ max 0\ max \ixiJt) - x^y^At) - y)\‘^ ] 

^,3=B,...,J I (x.y)enf.(t)" 


and finally 


^^sit) = - f yc^it,y)utit,y)dy + f utit,y)dy, 

dt Jniit) 

= miit)-c^it, 0 )ntgit) + o( max lyM . 

\y&ni{t) J 


(23) 
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To obtain the system of equations for dynamics of populations of couples we need 
to apply a similar treatment for masses and their locations. 




d rljw 


u^(t, X, y)dydx 


r^TW ^ r 

dt Jd 


lUt) 


lUiW 

lUt) 


u^{t,x,y)dy 


dx 


diT-iW 


ij-i(t) 


■ rd it) 
dii_iit) 


u%tJTit),y)dy 


u‘^{t,l'l^_^{t),y)dy 


ri^t) 

IdiW 


dtu‘'{t,x,y)dy 


+ ^ljit)u‘'{t,xjj(t)) - ^lj_^(t)u%t,x,l^j_^{t))jdx 

fit (t) fi^j (t) 

+ / u‘'{t,lTit),y)dy - / u‘'(t,l^^{t),y)dy 

= / dtu‘^{t,x,y)dx + / u^(t,xjUt))dx— / u^{t,x,l^_i{t))dx 

Jnf-it) JiiTW 

+ [ u‘'{tJT(.t)^y)dy- [ u‘^{t,l]^^{t),y)dy 

d^t it) Jo.!, it) 

= / dtu^{t,x,y)dx + / dxU^{t,x,y)dy + / dyu‘^{t,x,y)dx 

Jn-.it) Jn-.it) JQ-.it) 

= / [T{t,x,y) - c^{t,x,y)u’'{t,x,y)]dxdy. 

dn-^it) 

As previously, to differentiate the locations function {xiij,y^j){t) = ^^■^{x,y)u'^{t, x,y)dxdy 

we start with its first moment that is {x'ij,yij){t) = ^^s^{x,y)u'^{t,x,y)dxdy 




d r'-^d) pit it) 
dt Jifi^it) dif it) 


{x,y)u‘'{t,x,y)dydx 


/■C(‘) d 

fT w 

/ - 

/ ix,y)u‘'it,x,y)dy 


yT-id) 


dx 


dt 


IT it) 


Tit) 


T-tit) 


ilTit),yK{tjrit),y)dy 




it) 


'iTiW 


iiTTt),y)uTt,iT-iit),y)dy 
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+ 


+ 

+ 


/ {x,y)dtu‘'(t,x,y)dxdy+ {x,lUt))u‘'{t,x,lUt))dx 

- [ {x,l^j_^{t))u^t,xJ^^_^{t))dx 

JnY'd) 

[ iir{t),y)u%t,irit),y)dy- [ ^ ilT.,{t),yK{tJT-At),y)dy 
JnUt) JnUt) 


f 

/nj(t) 

/ {x,y)dtu‘'{t,x,y)dxdy 

/ {l,0)u‘'{t,x,y)dxdy + {x,y)dj;U^{t, x,y)dxdy 

/ {Q,l)u‘'{t,x,y)dxdy + / {x,y)dyu‘'{t,x,y)dxdy 

JnA (t) JnA (t) 

/ (a;, 2 /)[T(<,a:, y)-c‘'(t,a;, 2 /)M^(t,a;,y)](ixdy+ / ( 1 , x, y)(ia;dy. 

JnA(t) Jn9At) 


As the derivation of the closed form for couples goes in line with reasoning for the 
male and female populations, we apply equations (14) and (15). We omit details 
concerning approximation of particular functions, because similar calculations were 
performed in the case of male population. 

For a generic marriage function T(t,a;,y) we are not able to obtain the closed 
form, and most we can conclude is: 




-c‘=(t,x)j(t),y“ (t))m“j(t) + / T{t,x,y)dxdy 

J^Ut) 

{x,y)T{t,x,y)dxdy 


- iXii{t),yfAy)y{t,x‘r.{t),yfAt))m^At) + (1, l)mlAt) 


Hence, in further investigations we consider a specific marriage function given in 
equation (2). To obtain the desired form of ODEs describing dynamics of couples, 
we need to approximate the integral of the underlying function over cohort (t ): 


/ F{t,u'^{t,x),ud (t,x),u^{t,x,y))dxdy = (24) 

In<= (t) y)Hx)g(y) [■w"‘(i, x) - /“ u‘^(t, x, y)dy] [u-f(t, y) - x, y)dx] } dxdy 

7 + ir x) - u^(t, X, y)dy] dx + g(y) [uf{t, y) - /“ x, y)dx] dy 

To obtain the approximation of the first order with respect to the structure 
variables, we apply formulas (12)-(15). Observe that the denominator of function 
F does not depend on structure variables and hence, there is no need of integrating 
it over H“j (t). Let us start with the integral of the numerator. It consists of four 
parts which can be approximated separately: 
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Iniit) ^{x,y)h{x)g{y)vJ^ {t,x)u^ {t,y)dxdy = 


e{xT{t), y^^ {t))h{xT{t))g{yl {t))mT{t)mf (t) 

"'"^■‘'''2' ' '''^( max It/— r'^''^'''2 
y&nUt) 


+0(^max^) |x - a;™(t)|2) + 0{ max ^ \y - yj (t)|2) 


^OO 

^^^^y)Kx)g{y)u'^{t,x) J u''{t,z,y)dzdxdy = 

J 

^(^T{t):yvjit))HxT{t))g{yyj{t))mY'{t)mlj{t) + 


v^B 


+ max ol max | (x^j (t) - x, y^j (0 - 2/) N 

v=B,...,j y(x,y)ens^(t) j 

^ix,y)h{x)g{y)u^ {t,y) j u‘'(t,x, z)dzdxdy = 

J 

Vj it))HXiyjit))9{yj {i)Wj + 


u^B 


+ max O max |(x°^(i)-x,y°„(t)-2/)|2 

poo poo 

faf.(t) ^ix,y)h{x)g{y) J u^{t,x,z)dzj u^{t, z,y)dzdxdy = 

J 

v,'w—B 


+ max 0[ max |(xS^(i)-x,yS^(i)-y)|2 


To give the idea how the approximations are obtained, we present calculations 
of the third part: 

/ &ix,y)h{x)g{y)u^{t,y) u''{t,x, z)dzdxdy = 

JniAt) Jo 


iTd) 

iT-id) 


2X = 


iTit) 


noo nL- 

h{x) / u‘'{t,x,z)dz / Q{x,y)g{y)u^ {t,y)dyd: 

Jo Ji^At) 

pOO 

Kx) u%t,x,z)dz 0(x,y/(t))5(y/(t))mJ(t) + ||u^(t,-)|lLi(nJ(t)) •C’(|y-2//Wn 

(h{x) 0 (a:,y/(t))g( 2 // {t)Wj(t) + ||u^(i, •)llLi(nJ(t))C’(l2/“ 2//WH ) 


^T-id) 

^ rCCO riLit) 

J 

w—B 




a^B 


, z)dzdx 

-vA) 
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( max \{x'i^(t) - x,yl 

w—B 

^(2;LW)[0(a;L(i):2// {t))9{yj it)Wj W]mLW+ max O f max |(a;-^(t) - 

W—B ^ 

Let us define 

J 

v—B 

J 

+ 0(a;LW>2// {t))Kxtw{t))9i.yi{t)Wj (25) 

W—B 

J 

+ X! ®«™(^)>2/L(^))^(^ii«(^))5'(yL(^))"^L(^)”^L(i) 

v,w—B 

and observe that Nij(t) denotes the approximation of the numerator of formula 
(24). We can easily conclude that Nij{t), defined below, denotes an approximation 
the numerator of formula (24) multiplied by vector {x,y): 

Nij{t) = {x^{t),yUt))e{x^{t),yUt))h{x'^{t))g{yUt))m^{t)mUt) (26) 


V—B 

J 

+ Y (^iw(.t)^yj (.t)){^iXiwit),yj {t))Hxiwit))9{yj (0)mj(i)m-^(t) 


W—B 

J 


+ Y (®L(0:yL(^))®(®L(i),2/L(*))^(^L(i))5(ySj(0)mL(0mL(0- 

v,w—B 

The denominator of formula (24) is made up of two parts, which we approximate 
separately: 

poo r poo 1 •J J 

/ h{x) u"^{t,x)- u‘'{t,x,z)dz da; = ^/i(a;™(t))m7*(t)- ^ /i«^(t))m°j(t) 

^ 0 . ♦-'0 _ „• D „ X D 


i^B 


^,3—B 


i=B,...,J 


:^r2^ (t) 


+ max Oi max |a:™(t) — a:M + max Oi max |(a:° (f) — x, yL(t) — y)| 


9{y) 


i,j=B,...,j y(x,y)e%(t) ■' 

J J 


pOO 

(t, y)- u'=(t, 2 , y)dz dy=Y didj (0)mj(t)+ Y 9{ytiit))m%it) 

“'0 -I j=B z,j=B 


+ max oi max \y{(t)-yf]+ max o{ max \(xiJt) - x,y^Jt) - y)\ 
Let us write 

,/ .7 

Dij{t) = 7 + ^ h(x™(t))m™(t) - Y Hxij{t))m^ij{t)+ 

i—B i,j—B 


,{t)-y)\ 

x,yL{t) 


(27) 
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,/ J 

+ diVj it)Wj (t) - 

i—B i,j—B 

and understand that Dij(t) is an approximation of the denominator of formula (24). 
Now F(t,u^{t,x),u-f (t,x),u'^{t,x,y))dxdy is approximated by Y -[*) ; while 

/n“ is approximated by Y -[tj ■ 

Having the approximations, we obtain equations for mass and location of couples: 

s’"***' = 

+ 


-c%t,xlJt),y^Jt))m'=Jt) + 


Dij (t) 


+ h.o.t 


[( 1 : 1 ) - {x^ijit),ytj{y)V{t,x^,j{t),yfj{t))] 


Dij (t) 


+ h.o.t, 


(28) 


(29) 


where h.o.t stands for higher order terms, which are of second order. 

After neglecting higher order terms in equations (16)-(23), we are able to present 
the numerical scheme for masses and locations for the males, females and couples: 



= -c^{t,x^{t))mr{t), i = B + l,...,J, 


= 1, z = -S + l,...,</, 


= -c™(t,0)m^(t)-5,c™(t,0)n-(t) 


+ T.ij=B Z?"* (^: ^ij (1): yfj it))mij (t) , 


= m^(t)-c™(t,0)n^(t), 

x^(t) 

J 0ifTO^(<)=0, 

1 rh|(t)t otherwise, 


= -cl(t, ijj {t), j = B + 


= 1, j = 5-1-1,..., J, 


= -cl(t, 0)rh^g(t) - dycf{t, O)n^(t) 

tMH) 

+ T,ij=B W’ 

= miit) - 

yiit) 

J 0 if = 0, 

1 


(30) 


(31) 


(*Ei,(i),i'y(i)) 


-c%t,x'r^{t),yfj{t))m^,j{t) + i,j = 

[( 1 : 1 ) - ix^ijit),yijiyW{t,Xijit),y^jit))] rh’i^it) + 

f 0ifmUt)=0, 

i otherwise, 


Bij (t) 


(32) 
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where Nij{t) and Dij{t) are given by (25), (26) and (27), respectively, except 

that all the locations of the cohorts and their masses are replaced with relevant 
variables with hats. 

Observe that the right hand side of the last but one equation, e.i. for ^ (t ), yfj (t)), 

consists of terms only of the form 

Above differential equations need to be supplemented with the initial condi¬ 
tions at the beginning of the interval [t„,t„+i). Let us start with ODEs for 
the males and females. For equations evolving in the interval [tQ,ti) the initial 
conditions for our scheme a;f“(0), to"*( 0), yj(0), mj(0), for i,j = are 

known based on the initial conditions of the original problem, while th^^(O) = 0 
and n^p(O) = 0. For the corresponding ODEs evolving in [tn,tn+i), where 0 < 
n < TV, we take xY"{tn) = limj_,.j-i™(t), m^{tn) = mY"{t), yj{tn) = 

yj(t), rh{{tn) = m{{t), for i,j = iT„,..., J and as previously 

= finjin) = 0 and = n^„(tn) = 0. 

In the case of couples, equations evolving in the interval [to,ti) are supple¬ 
mented with the following initial conditions: TOg^^^(O) = Tn‘[g^{Q) = = 0, 

= (O’O), = (0,0), for j = Bo,...,J and 

(0)> yfBo (0)) = (0, 0), for * = So,..., J. 


3. Well-posedness and convergence of the method. 


3.1. Well-posedness of the structured population model. Existence, unique¬ 
ness and Lipschitz continuity with respect to initial data and model parameters 
of a nonlinear version of model (1) were established in ref. [17] using the frame¬ 
work of nonnegative Radon measures (A1"'"(R+) x Ad+(R+) x Ad+(R^),(i), for 
d = di + di + d 2 with Lipschitz bounded distance 


di{yi, ^ 2 ) = sup < / ipd{fii - ^ 2 ) ■ G C (RT,_;R) and 




<U, i = l,2. 


It was shown in ref. [17] that, under the assumption that 

c^c™,/3/,^'",h,gGC([0,T];Wi'-(M+)) 

and 

c^0 G C([0,T];Wi’°°(R^)), 

a Lipschitz semigroup of solutions is generated in the sense of Definition 1.1. 


3.2. Solvability of the EBT algorithm. Proof of solvability of the EBT method, 
i.e. existence and uniqueness of solutions of system (30)-(31), is not straightforward, 
because of the specihc definition of the dynamics of the cohorts involving quotients 
of the variables with denominator which is not necessarily separated from zero. 
Similar problem for the EBT algorithm for a scalar structured population model 
was solved in ref. [6]. The strategy proposed in ref. [6] can be applied also in our 
case. 

It consists in considering the problem in a restricted domain, in which the right- 
hand side is locally bounded and locally Lipschitz-continuous, which yields local 
in time existence and uniqueness of solutions. The main difficulty is in showing 
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that the solutions remain in the restricted domain locally in time, i.e. the values of 
yfjit)) remain in a cone 

(0,0) < < {Ci{t),C2{t)) (33) 

for certain Ci(t), i = 1, 2, which is equivalent to proving the following inequalities: 


0 < 



cyt). 


(34) 

(35) 


While the first inequality is obvious, the second one needs justification. Let us 

X^(t) 

remind that x^j{t) = and take Ci{t) = t + C, where C is a constant. Then, 

the inequality (35) is equivalent to {t) —t < C. 

Equation (28) yields 


1 


Taking 


we obtain 


4(t)-^ j -xiymyt) 

{xy{t) - x^,j{t))Q{xy{t),yf{t))h{xy{t))g{y^{t))my {t)rhf {t) 

J 

+ {xy{t) - x'lyt)) &{xy{t),yy{t))h{xy{t))g{yy{t))my{t)mlj{t) 

v—B 

J 

w—B 

J 

+ Y (^L(i) -^o(O)0(^^L(O:ySj(O)^(^L(O)5(ySj(O)™Sj(O™L 

v,w—B 

^imax(0 = argmax^^B jx)=^(t), 


1 


^ (^?max(i) -t)<0 for > xy{t) =C + t. 

Note that ^ (x^i^^bt) - t) = ^ (^?max(0 7 ^r(O) and therefore \x^i^^yt)-xy{t)\+ < 
l^imax(O) ~ ir(0)l'''- Finaly 0 < a;°j(t) < C + t. For further details of the proof, we 
refer to ref. [ 6 ]. 


3.3. Remarks on the convergence of the method. The essential step in the 
proof of the convergence of the numerical method is the stability result for the Lips- 
chitz semigroup presented in ref. [17], where the metric space {E, p) from Definition 
1.1 is understood as (A4+(R-|-) x Af'*'(R-|-) x Af+(K+),(i), for d = di +di +^2 and 

[ ipd{pi - P 2 ) ■ & C^(R!,_;R) and ||i^||wi,~ < 1 

By 

As reported above, equation ( 1 ) is associated with a Lipschitz semigroup, [17]. 
Therefore, due to Proposition 1, the remaining step in the proof of convergence 



di{pi,P2) = sup' 
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boils down to the estimate of the tangential condition, which is the right hand side 
of the estimate (8). To identify the missing part of the proof, let us define two 
continuous operators, projection and extension, in the following way: 

P : {M + (R+) X M + (R+) X Ad + (R^), d) M4(J-Bo+n+l) ^ ^3(J-Bo+n+lf 

E : k4(J-Bo+"+i) X (X+(R^) x Af+(R+) X Af+(R2 


given by the formula 

Ei{xT{t),mT{t)), {y( (t)), 

J J J .1 

i=B{t) j=B{t) i=B(t) j=B(t) 

One observes that P o £1 = Id, while the E o P does not need to be the identity. 

From the derivation of the method we obtain immediately an estimate for the 
tangential condition between P(u™(t, •), •), •, •)) and 

{xY'{t),mY"{t)), {y( (t),mf(t)), {x‘ij{t),yfj{t),m^j{t)) in the finite dimensional space 
of the approximation, with the appropriate accuracy. To present the rigorous proof 
of the convergence of the method we need the estimate of the condition of tangen- 
tiality between (u™(t, •), w-f (t, •), •, •)) and 

iJ2i=B{t) '’^Ti^)hxY'{t)},J2j=B{t) (4)}I Z]i=B(t) J2j=B{t) 

in the space of measures equipped with an appropriate distance. Similar proof was 
carried out in ref. [ 6 ]. Due to the lengthy of technical calculations needed in this 
step, a complete proof is deferred to a forthcoming paper. 
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